Astronomy & Astrophysics manuscript no. aal0810 


© ESO 2009 


January 27, 2009 





ON 
O 

o 

(N 
G 

(N 



Oh 
6 

C/3 



> 

(N 
(N 

o 

On 
O 



13 



Notes on the disentangling of spectra 

I. Enhancement in precision 

P. Hadrava 



Astronomical Institute, Academy of Sciences, Bocm II 1401, CZ - 141 31 Praha 4, Czech Republic 
e-mail: had@sunstel . asu .cas.cz 



Received 15 August 2008 / Accepted 16 October 2008 



ABSTRACT 



Context. The technique of disentangling has been applied to numerous high-precision studies of spectroscopic binaries and multiple 
stars. Although, its possibilities have not yet been fully understood and exploited. 

Aims. Theoretical background aspects of the method, its latest improvements and hints for its use in practice are explained in this 
series of papers. 

Methods. In this first paper of the series, we discuss spectral-resolution limitations due to a discrete representation of the observed 
spectra and introduce a new method how to achieve a precision higher than the step of input-data binning. 

Results. Based on this principle, the latest version of the KOREL code for Fourier disentangling achieves an increase in precision for 
an order of magnitude. 
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c3 1. Introduction 



Disentangling spectra of binary and multiple stars enables us to 
determine efficiently the orbital parameters and simultaneously 
to separate the spectra of the component stars. This numerical 
technique performed in either the wavelength domain (Simon 
and Sturm 1994) or its Fourier image (Hadrava 1995) has been 
applied successfully in numerous studies of individual stellar 
systems. However, some users failed in their attempts or were 
unable to take full advantage of the method, partly due to a mis- 
understanding of its principles. A review of the Fourier disen- 
tangling has been provided by Hadrava (2004) together with the 
release in 2004 of the author's code KOREL, but, regarding new 
improvements of the method, this review is already out of date. 
The purpose of the present series of papers is to explain some 
common mistakes, provide practical hints for using the method, 
and present its new developments. 

In this paper, consequences of the discretization of the ob- 
served spectra are discussed in Sect. [2] A new method for en- 
hancement of the spectral resolution in disentangling spectra is 
introduced in Sect. [3] Results and their implications are briefly 
summarized in Sect. [4] 



2. Sampling of the input spectra 

In their study of disentangling, Hensberge et al. (2008) specu- 
lated about "expense" at which the computational efficiency of 
the spectral method dealing with the Fourier image surpasses 
the method of singular-value decomposition in the wavelength 
domain. They suggested that, among other things, it may be the 
need of having the input observed spectra sampled on a common 
grid equidistant in the logarithmic wavelength scale 



where Aq is an arbitrarily chosen reference wavelength. It should 
be noted that the same assumption is commonly imposed on the 
solution in the wavelength domain as well, as it has been de- 
scribed by Simon and Sturm (1994) and it is also obvious from 
the explanation in Fig. 1 of Hensberge et al. (2008) or their ex- 
ample in Appendix of the same paper. 

The uniform sampling of input data simplifies the solution, 
but this assumption may be avoided in both methods of disen- 
tangling. In the Fourier view, it is obvious that the Fourier trans- 
forms 7(y) of the observed spectra I(x) in chosen (equidistant) 
sampling frequencies may be calculated directly according to 
the definition 



f I(x 



) exp(iyiiX)dx 



(2) 



from any original (even non-equidistant) binning xi if the func- 
tion I(x) is suitably interpolated, e.g. by the simple linear for- 
mula 



I(x) 



X[-\-\ — X X — X[ 



xi+i - xi 



xi+\ - x. 



(3) 



for x € (xi,xi +i ) 



x — c IuA/Aq 



(1) 



The common practice of interpolating I(x) first to the equidis- 
tant grid points and then using the Fast Fourier Transform saves 
the computer time (at some expense of accuracy), but is not in- 
evitable in Fourier disentangling. 

In the wavelength-domain solution, the single off-diagonal 
matrices (N in notation of Simon and Sturm 1994) shifting the 
spectra of component stars to their appropriate positions in indi- 
vidual exposures may be replaced by wider band matrices if the 
observed spectra are not sampled in the same equidistant set of 
the logarithmic wavelengths into which the component spectra 
are to be separated (Simon and Sturm 1994, p. 287: "The sub- 
matrices of M, Nai and Ng,-, are rectangular band matrices with 
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a bandwidth depending on the differences in dispersion of the 
wavelength scales of c and x"). A simple possibility is to use a 
matrix N with two non-zero elements in each column given by 
Eq. (0. This is, however, again equivalent to a suitable resam- 
pling, and a subsequent solution in the convenient representa- 
tion of the observed spectra. An additional significance of using 
bandwith matrices to refine the model is discussed in the next 
Section. 

Although any resampling implies some smoothing of the in- 
put signal and thus a loss of information in the high-frequency 
modes, it is inevitable provided the observations are not directly 
obtained in the required data bins. The question is therefore not 
about performing a resampling, but how it can be performed 
best. This problem, which is common to any method of disen- 
tangling, is related to the more general task of optimal data- 
processing of observed spectra (cf. Hensberge 2004) and its as- 
pects in disentangling will be studied in detail elsewhere. 

A consequence of discretizing observed spectra is limitation 
of the accuracy at which the radial velocities are determined. 
Until now, a common practice in both the wavelength-domain 
and Fourier disentangling was to round the expected Doppler 
shift to an integer multiple of the radial-velocity step. This lim- 
ited resolution of Doppler shifts in the individual spectra also 
limits the precision of the disentangled orbital elements and 
the sharpness of the separated spectra of component stars. A 
straightforward means of improvement appears to be a choice 
of a smaller sampling step. However, the resolution is limited 
by the detectors in any case. On the other hand, it is intuitively 
evident that if we have a set of spectra with mutually shifted 
sampling, we can also reconstruct details on a sub-pixel scale. 
For instance, if a very narrow line with a sub-pixel width moves 
from a given pixel in some exposures to the neighbouring pixel 
in other exposures, then its position may be found precisely from 
the time and its width from the duration of this transition. It is 
thus worth investigating the limits of resolution in details. 

3. An increase of spectral resolution 

In my spectral method, the shift in the spectra I(x) of each com- 
ponent (which we wish to separate from the observed superposi- 
tions), in the logarithmic wavelength scale x defined by Eq. ([1) 
for a value v of the radial velocity, is given by the convolution of 
the spectra with the shifted Dirac delta-function S(x - v), 



/'(*) = I(x - v) = /(*) * S(x - v) , 

which implies, in the Fourier transform (x 
by a function exp(iyv), 

I'iy) = I(y) exp(ryv) 



(4) 

y), a multiplication 



(5) 



(cf. Ij in Eqs. (1) and (2) of Hadrava 1995). This simple expo- 
nential function can be evaluated precisely at each frequency y. 
However, due to the limited number N of the modes taken into 
account, its inverse Fourier transform will generally produce a 
wider peak with some ghosts on its sides resembling interference 
fringes. Only in the special case of v being an integer multiple 
of the grid step, the period of function exp(iyv) is in resonance 
with the interval length in ^-representation and a sharp shifted 5- 
function coinciding with a grid point of the x-representation can 
be reproduced. For that reason the radial velocity was rounded 
to the nearest grid point in the Fourier disentangling also and it 
explains why the radial velocities or their residuals calculated 
by the original KOREL-code were quantized depending on the 
radial-velocity step. 




Fig. 1. Discretization of a Lorentzian profile I(x) (the smooth 
thick line) centered on the pixel position should yield a symmet- 
ric distribution of counts D[I](x) in neighbouring bins (the thick 
step function). A slightly shifted profile I'(x) (for 0.2 pixel- width 
in this figure - see the thin lines) results in an asymmetry of the 
counts D[I'](x), which in turn enables us to determine the line 
position at a precision below the pixel width. 



However, owing to the resolution in the digitalized values 
of intensity read from individual detector pixels, the position of 
spectral lines wider than the sampling step can be deduced with 
an accuracy exceeding the step width (cf. Fig[TJ. Alternatively to 
a convolution with the shifted ^-function, a shift of a spectrum 
I(x) for value v can be expressed as a Taylor expansion 



oo 

I(x - v) = ^ -I U \x)(-v) J . 

7=0 J ' 



(6) 



which usually converges rapidly for small values of v. In a dis- 
crete equidistant representation xi with the step A x , the first two 
derivatives can be approximated by finite differences 



1 



2A X 
1 



I (z >(x k ) - -r^(Kx k+1 ) - 2I(x k ) + /(**_!)) . 
At 



(7) 
(8) 



Therefore in the vicinity of the grid point Xk a small shift of / can 
be expressed in terms of values in this and the two neighbouring 
points as 



/fot-v) - I{x k )-I m {x k )v+ l -I {2 \x k )v 2 +o(v') 

V 

- I(xk) - j^-(Kxk+i) - I(x k -i)) + 
+ 2A 2 " (/(X * +l) ~ 2I<yXk) + I{Xk - l)) + ° {V ) ■ 



(9) 
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Fig. 2. Fourier transform of the discretized shifted profile D[I'] 
corresponds to the Fourier transform of the centered profile D[I] 
multiplied by a correction, the real and imaginary parts of which 
are drawn by the upper and lower (distorted sinusoidal) thick 
lines, resp., while their approximations (fTTT i are drawn by the 
thin sinusoidal lines. 



This implies that the operator S(x— v) of the shift is approximated 

by 

S(x - v) * S(x) - 2T-(<*(* + A *) - - A *» + ( 10 > 

v 2 , 
+ 2Af (<5(X + x) ~ 26(X) + 6(X ~ Ax)) + ° ( } 

and its Fourier transform 

v 

expO'yv) a 1 - _(exp(-ryA x ) - exp(;>A x )) + 
2A X 

v 2 , 
+ ^2(exp(-!>A A .) - 2 + expO>A x )) + o(v J ) = 

2 

= 1 + 1- sb(yA,) + ^(cosCyA,) - 1) + o(v 3 ) . (11) 
A x A^ 

It can be seen in Fig.|2]that a ratio of the two Fourier trans- 
forms of mutually shifted profiles is approximated well by this 
simple sinusoid for small values of y, while for the higher- 
frequency modes (drawn closer to the middle of the figure) 
higher harmonics contribute significantly. This is an obvious 
consequence of the fact that the approximations ((T) and (JHJ of 
the derivatives are more accurate for the lower modes, which do 
not change significantly on the scale of A*. The exact shape of 
the shift operator depends on the spectrum to be shifted, unless 
v is an integer multiple of A x , It means that the value of v cannot 
be disentangled with unlimited precision from roughly sampled, 
unknown spectra. However, already the application of the cor- 
rection ( fTTT i improves the precision of the disentangling signifi- 
cantly, and the accuracy could be even higher for disentangling 
constrained by a template spectrum. 



The above described procedure of reconstructing component 
spectra from a large set of observations should not be confused 
with a simple interpolation (given e.g. by Eq. OJ) in a single ob- 
servation or between grid points of some theoretical models. For 
instance, in model atmospheres, the dependent variable (e.g. the 
specific intensity) is calculated usually for chosen exact values of 
the independent variable (wavelength) from which it can be in- 
terpolated to other values or integrated over some regions of the 
independent variable. On the other hand, in true observations, 
the values read at individual detector pixels are the quantities 
integrated over some interval of the wavelength, which provide 
some constraint only on the inner distribution. Without any addi- 
tional information, these values may be used in a single exposure 
as an estimate of the variable for the middle of the interval, while 
closer to its edges the value interpolated between the neighbour- 
ing bins is more appropriate. However, for a set of exposures, 
the information can be combined to reveal partly also the sub- 
pixel structure, or, using the above described procedure, to find 
subpixel mutual shifts between the exposures. 

As an example, we show a result of disentangling of simu- 
lated data. Twenty spectra uniformly covering one period (which 
is taken to be a unit of time) of a double-line binary on circular 
orbit with chosen radial-velocity semi-amplitudes K\ = 50km/s 
for the primary, and Ki = lOOkm/s for the secondary were cal- 
culated. For each component, one line with a Lorentzian pro- 
file (with central depths 0.3 and 0.2 of the common continuum 
and half-widths equal to 30 and 40 km/s for the primary and 
secondary, respectively) was included. A pseudo-random noise 
scaled to amplitudes n = 0%, 0.5%, 1.%, or 2.% of the con- 
tinuum level was added. The spectra were sampled by integrat- 
ing in bins of width corresponding to 10 km/s. The results of 
disentangling obtained using the KOREL code in its old ver- 
sion (KOREL04 released by the author in 2004) and in its new 
version (KOREL08) of enhanced precision, are compared in 
Table [TJ In this Table, S denotes the integrated square of spec- 
tra residuals, and A7o is the difference in units of the period be- 
tween the solved epoch of periastron (defined by fixed periastron 
longitude) and its true value Tq — chosen for the simulation. 
Similarly, AKi are the differences between the calculated and 
true radial-velocity semi-amplitudes of the primary and Aq for 
the mass ratio (q = Mz/M\ = Ki/K 2 = 0.5). 

It can be seen from the results that the squares S of the resid- 
uals consist of a part approximately proportional to the square 
n 2 of the noise, as can be supposed, but also an other addi- 
tive, almost constant part, which is comparable to the 1 % noise 
in the solution with the classical KOREL04, but is suppressed 
for at least two orders in the super-resolution KOREL08. This 
part is obviously due to the discrepancies between sampling of 
the component spectra in different exposures shifted by a non- 
integer multiple of the sampling step. This contribution depends 
on the shape of the spectrum and its importance on the level 
of the noise. This explains why in preliminary applications to 
real data the new method yielded significantly superior results in 
some cases, but only a negligible improvement in other cases. 

Similarly, the errors in orbital parameters have a part that 
increases with the noise and a noise-independent part, which is 
significantly smaller in the solution based on the new KOREL08. 

4. Conclusions 

The correction provided in Eq. (fTTT i for the residual part of the 
radial velocities over an integer multiple of the sampling step 
improves the Fourier disentangling significantly. With this re- 
sult and other improvements completed by the author to recent 
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Table 1. Comparison of the old and new KOREL solutions. 



Parameter 


n[%] 


KOREL04 


KOREL08 


c 


u 


3 £ 1 


U.Uiz 




0.5 


4.98 


1.22 




L. 


8.76 


4.85 




2. 


23.59 


19.33 


AT 1 / ) ■ JT 

Alo [rerwa\ 


U 


4.y x to 


4. X 10 




0.5 


4.9 x 10~ 4 


6.6 x 10~ 5 




1 


4.9 x 10~ 4 


2.3 x 10 -4 










AK[ [km/s] 





0.49 


0.07 




0.5 


0.49 


0.03 




L. 


0.49 


0.06 




2. 


0.33 


0.18 


Aq 





0.010 


0.0005 




0.5 


0.010 


0.0002 




1. 


0.010 


0.0014 




2. 


0.013 


0.0040 



versions of the KOREL code, the version of 2004 is no longer 
supported, and we recommend using for true applications the 
version of 2008. 

Due to the equivalence of Eqs. (0]l and (Q, the wavelength- 
domain solution could be improved similarly if the single off- 
diagonal matrix N would be replaced by a three- (off-)diagonal 
matrix fllOt . or even by a more complicated matrix, if it should 
also include an interpolation from a non-uniform sampling of the 
input data. A possibility for using band-matrices was mentioned 
by Simon and Sturm (1994), and in more detail explained by 
Sturm (1994) (cf. also Hensberge et al., 2008). 

The method described here could be improved to achieve 
an even higher precision for known component spectra (i.e. for 
the constrained disentangling), for which the higher harmon- 
ics of the approximation given by Eq. ( TTTT i could be estimated. 
Analogous numerical refinement either in direct or Fourier space 
could be useful also in other methods in spectroscopy (e.g. in 
methods using the broadening function) as well as in data pro- 
cessing in other fields of astrophysics. 
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